\documentclass[12pt]{article}
\usepackage{bm}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{booktabs}
\usepackage{float}

\begin{document}
\section{Introduction}
Various fluid phenomenons are around us in daily life. From Ball games to automobiles, ships and airplanes. However, without good scientific understand of it, we can not predict the behaviour of them. And nowadays, it is important to have a deep insight into the understanding of fluid flow from microscopic scale (Oil, water and gas in reservoir) to macroscopic (Aerodynamics for air planes). Better understanding of fluid could help us to predict the behaviour in order to improve many kinds of technical devices, or make them work more efficiently. And normally, the fluid phenomena are coupled with the interactions with other physical, chemical and biological processes. These multi-physics flow problems are also very common in our daily life.\\

In general, we have three ways to understand and predict the behaviour of fluid flow: experiment, analytical solution of the underlying equations and numerical simulation.Due to the complicity of fluid equations, analytical solution are archived for only very few fluid flow problems. Experiments help us to understand the fundamental principles of fluid, explore the complex flow. However, these experiments are extremely time-consuming, and highly specific equipment which are very expensive are required. And for some particular problems like weather prediction, it's not possible to implement an experiment for it. Therefor, numerical simulations becomes important for fluid dynamic research.With the rapidly development of computer power and efficient numerical methods, numerical simulations are capable for more and more fluid phenomenas, and help us better understand the fluid and improve our design.\\

In this paper, we are going to introduce a novel method for fluid simulation Lattice Boltzmann Method, which discrete the Boltzmann equations and solve it to simulate the fluid flow instead of solving the Navier–Stokes equations. It has several advantages over other CFD technique especially in simulation of flow with complicated boundaries, multi-component multi-phase flow simulation, evaporation, condensation, cavitation as well as high Reynolds number flows.\\

\section{Overview of Lattice Boltzmann Method}
\subsection{History and introduction}
Lattice Boltzmann Method originated from Lattice gas automata (LBA). LBA models fictitious particles dynamics with discrete space, velocities and time in a simplified way. All the gas particles are located on the lattices and  with certain velocities.At each time step, two procedures are carried out, streaming and collision. In streaming, each particle will move to the neighbour lattices site according to the velocity it has. In collision step, the particles velocity will be determined from itself and the neighbour particles with collision rules. The state of each lattice site is Boolean. This method has several major disadvantages such as lack of lack of Galilean invariance, statistical noise which limit its usage.\\


In order to reduce the statistical noise from Boolean calculation, density distribution function was introduced to replace the Boolean variable. Due to the replacement of Boolean, the collision rule is also replaced by continuous collision operator.\\

Instead of solving the conservation equations of macroscopic properties (mass, momentum, and energy) numerically, LBM models the dynamics of fictive particles, and such particles perform consecutive propagation and collision processes over a regluar lattice mesh. Then fluid motion could be obtained without solving the partial differential equations.\\

\subsection{Kinetic theory and Boltzmann equation}
We simplify the gas into hard spherical particles with same size and mass and  there are only elastic collisions between each other.With the information of position and velocity of each particles, together with classic mechanical theories, it would be possible to predict the future state of the system.\\

A function $f(x,\xi,t)$ is introduced to describe the distribution of particles. It represents the number density of particles at the time $t$, and in the position x with the velocity $\xi$. We consider the gas in a control volume $dV=[x, x+dx]$. And an external force $F=ma$ is applied to the system. At time $t$, and the number of particles with the velocity $[\xi, \xi+d\xi]$ is $dN=f(x,\xi,t)d\xi dV$. After the time $dt$, if there are not collision, the location of these particles will be $x'=x+\xi dt$, and the velocity will be $\xi'=\xi+adt$. Therefore, we have equation:

\begin{eqnarray}
f(x+\xi dt,\xi+d\xi,t+dt)d\xi dxV-f(x,\xi,t)d\xi dx=0
\end{eqnarray} 

Or

\begin{equation}
\frac{\partial f}{\partial t}+\xi \cdot \nabla_{x}f+ a \cdot \nabla_{\xi}f=0
\end{equation}

Collision could result in change of the velocity of particles, so the number of particles in the control volume $[x, x+dx] \times [\xi, \xi+d\xi]$ and $[x', x'+dx] \times [\xi', \xi'+d\xi]$ will not be same any more. A collision term 
$\Omega (f)$ is introduced to account for the distribution change of particles. For ideal gas, the collision term could be described as:

\begin{equation}
\Omega(f(\xi))=\int [f' f_1' -ff_1]B(\theta, |V|)d \theta d\epsilon d\xi_1
\end{equation} 

$f(x,\xi,t)$ is the number density before the collision, and $f'(x,\xi',t)$ is the density after collision. $V=\xi_1-\xi$ is the relative velocity between two particles, $B(\theta, |V|)$ is a non negative function related with interactions between particles.\\

Function $\psi(\xi)$ is an arbitrarily function about $\xi$. If we make $\psi=1,\xi,|\xi^2|$, we can proof that:

\begin{equation}
\int \Omega(f)\psi(\xi)d \xi =0
\label{collision invariant}
\end{equation}

Any function which satisfy the equation \ref{collision invariant} is called Collision invariant. 

\subsection{H-theorem}
In 1872, H-theorem was introduced by Boltzmann to describes the increase in the entropy of an ideal gas in an irreversible process. A function called H function was defined as:

\begin{equation}
H(t)=\overline{lnf}=\frac{\int flnf d\xi}{\int fd\xi}=\frac{1}{n}\int flnf d\xi
\end{equation}

He proofed that The H function is a monotone decreasing function along with time:
\begin{equation}
\frac{\partial H}{\partial t} \le 0
\end{equation}
When it reach a minimized value, the system reach equilibrium:
\begin{equation}
 \frac{\partial H}{\partial t} =0
\end{equation}


\subsection{Maxwell distribution}
In 1860, Maxwell gave a probability equilibrium distribution of speed for particles in gas. It gives us the probability of a particle's speed being near a specified value in the function of temperature, velocity, and mass.

\begin{equation}
f=n\frac{1}{(2\pi R_gT)^{2/3}}exp[-\frac{(\xi-u)^2}{2R_gT}]
\end{equation}

Where $R_g$ is gas constant, $T$ is the thermodynamic temperature.\\


\subsection{Boltzmann-BGK equation}
In 1954, Bhatnagar, Gross and Krook \cite{BGK1954} introduced an approximation term $\Omega_f$ for the collision term $\Omega(f)$. And they proof that a simplified term $\Omega_f$ which replace the collision term $\Omega(f)$ should satisfy the following two properties:

\begin{itemize}
\item For collision invariant $\psi=(m,m\xi,\frac{1}{2}m\xi^2)$, equation
      	\begin{equation}
	\int \psi \Omega_f d\xi=0
	\label{BGK1}
	\end{equation}
       should be satisfied.
\item And it should satisfy Boltzmann H thereom:
	\begin{equation}
	\int(1+lnf)\Omega_fd\xi \le 0
	\label{BGK2}
	\end{equation}
\end{itemize}

They obtained the simplified term with the idea: The collision will lead the system to equilibrium distribution $f^{eq}$. And the change rate is proportional to the difference of $f^{eq}$ and $f$. And the scale factor $\nu$ is a constant. The approximation term $\Omega_f$ is obtained:

\begin{equation} 
\Omega_f=\nu[f^{eq}(x,\xi)-f(x,\xi,t)]
\end{equation}

The Boltzmann equation is simplified to:
\begin{equation}
\frac{\partial f}{\partial t}+\xi \cdot \nabla_{x}f+ a \cdot \nabla_{\xi}f=\nu (f^{eq}-f)
\label{Boltzmann-BGK}
\end{equation}

Equation \ref{Boltzmann-BGK} is called Boltzmann-BGK equation. And Bhatnagar, Gross and Krook also proof $\Omega_f$ satisfy the property \ref{BGK1} and \ref{BGK2}. A collision time $\tau=\frac{1}{\nu}$ also called relaxation time indicating the time interval between two collision is introduced in the equation.And Boltzmann-BGK equation is given as:

\begin{equation}
\frac{\partial f}{\partial t}+\xi \cdot \nabla_{x}f+ a \cdot \nabla_{\xi}f=-\frac{1}{\tau}(f-f^{eq})
\end{equation}

The macroscopic fluid density, velocity and energy could be calculated from the microscopic distribution function:
\begin{eqnarray}
&\rho(x,t)=mn(x,t)=\int f(x,\xi,t)d\xi & \\
&n\mathbf{u}(x,t)=\int \xi (x,\xi,t) d\xi&\\
&nR_gT(x,t)=\frac{1}{D} \int (v-\mathbf{u})^2f(x,\xi,t)d\xi&
\end{eqnarray}

Where $m$ is the particle mass, $D$ is the dimension of the space, $\rho,T$ is the macroscopic density and temperature. By applying Chapman-Enskog expansion, the macroscopic equations for mass, momentum and energy can be derived from the Boltzmann equations.The bulk viscosity is derived as:
\begin{equation}
\nu=\frac{\tau R_g T}{m}
\end{equation}

And the equation of state (EOS) relating pressure and density is given by:
\begin{equation}
p=nR_gT
\end{equation}


\subsection{Lattice Boltzmann Method}
Lattice Boltzmann Method is a special discretization of Boltzmann-BGK equation. Descretization of space, velocity and time are carried out in LBM. And this procedure greatly simplified the original Boltzmann equation. The location of particles in the space is restrained on the nodes of lattice grid, and the particle velocity is simplified into a very limited number of lattice velocities. We take a 2-D model as an example. This model is well known and widely used in the application. It is two dimensional and contain 9 velocities with the name D2Q9. This model is proposed by Qian et al\cite{qian1992}. In LBM, we assume all the particles has the same uniform mass (normally 1 mass unit is taken for simplicity). And the lattice unit($lu$) and time steps ($ts$) is important length and time unit in LBM. We only discuss with uniform mesh in this chapter ($\Delta x= \Delta y$).


\begin{figure}[!hpg]
\begin{center}
\input{d2q9.pstex_t}
\end{center}
\caption{D2Q9 lattice and velocities}\label{d2q9}
\end{figure}

The figure \ref{d2q9} shows the discretized velocity space $\{ \mathbf{e}_i \} ,\quad (i=0..8)$. The lattice velocity could be written as:

\begin{equation}
\mathbf{e}=e\begin{bmatrix}
0 & 1 & 0 &-1 &0 &1& -1& -1& 1 \\
0 &0 &1 &0 &-1& 1& 1& -1& -1 
\end{bmatrix}
\end{equation}

Where $e=\Delta x / \Delta y$ is the local lattice speed, and have the relation with sound speed as: $c_s=\frac{c}{\sqrt{3}}$
Similarly, D3Q19 model used for 3D LBM (Fig.\ref{d3q19}) has 19 lattice velocities and they are:

\begin{equation}
\setcounter{MaxMatrixCols}{20}
\mathbf{c}=c \begin{bmatrix}
0 & 1& -1& 0& 0& 0& 0& 1& -1& 0& -1& 0& -1& -1& 1 & 0& 0& 0& 0 \\
0 &0 &0 &1 & -1& 0 &0 &1 &-1 & -1& 1& 0 &0 &0 & 0& 1 &-1 &0 &-1 \\
0 &0 & 0& 0 & 0& 1& -1 &0 &0 & 0& 0& 1 &-1 & 0& -1& 0& -1& -1& 1
\end{bmatrix}
\end{equation}
\setcounter{MaxMatrixCols}{10}


\begin{figure}[!hpg]
\begin{center}
\includegraphics[width=2in]{d3q19.eps}
\end{center}
\caption{D3Q19 lattice and velocities}\label{d3q19}
\end{figure}



Correspondingly, the continuous distribution function associated with velocity is written as $f_i(x,t),\quad (i=0..8)$. We can obtain the Lattice Boltzmann Equation with D2Q9 model (Single Relaxation Time BGK):

\begin{equation}
f_i(\mathbf{x}+\mathbf{e}_i \Delta t,t+\Delta t)=f_i(\mathbf{x},t)-\frac{f_i(\mathbf{x},t)-f^{eq}_i(\mathbf{x},t)}{\tau}
\end{equation} 

As mentioned in previous chapter, colllision of the particles could be considered as a relaxation process towards a equilibrium. The equilibrium distribution function is defined as:

\begin{equation}
f^{eq}_i(x)=w_i\rho (x)[1+3\frac{\mathbf{e}_i \cdot \mathbf{u}}{e^2}+\frac{9(\mathbf{e}_i \cdot \mathbf{u})^2}{2e^4}-\frac{3\mathbf{u}^2}{c^2}]
\end{equation}

Where the weight coefficients for D2Q9 model are:
\begin{equation}
w_i=\begin{cases}
4/9 & i=0 \\
1/9 & i=1..4 \\
1/36 & i=5..8
\end{cases}
\end{equation}

For D3Q19 model, weight coefficients are:
\begin{equation}
w_i=\begin{cases}
1/3 & i=0 \\
1/18 & i=1..6 \\
1/36 & i=7..18
\end{cases}
\end{equation}

Macroscopic transport equations for mass, momentum and energy can be derived from Boltzmann equations using Chapman-Enskog expansion.\cite{ekexpansion} Kinematic viscosity $\nu$ in the D2Q9 model is obtained as:
\begin{equation}
\nu=c^2_s(\tau-\frac{1}{2})\Delta t
\end{equation}

Its unit are $lu^2ts^{-1}$. Note that $\tau>1/2$  for positive viscosity. Numerical difficulties can arise as $\tau$ approaches $1/2$. A value of $\tau=1$ is safest and leads to $\nu=1/6 lu^2ts^{-1}$.\cite{LBMODELING} The pressure is given by equation of state for ideal gas:
\begin{equation}
p=\rho c^2_s
\end{equation}


To implement a Lattice Boltzmann simulation, four major steps should be included in the code:
\begin{itemize}
\item Initialisation of distribution function $f_i(\mathbf{x},0)$
\item Collision step
	\begin{equation}
	f'_i(\mathbf{x},t)=f_i(\mathbf{x},t)+\Omega_i(\mathbf{x},t)
	\end{equation}
\item Steaming step
	\begin{equation}
	f_i(\mathbf{x}+\mathbf{c}_i\Delta t,t+\Delta t)=f'_i(\mathbf{x},t)
	\end{equation}
\item Macroscopic hydrodynamic quantities computing
	\begin{eqnarray}
	&\rho(\mathbf{x},t)=\sum_{i} f_i(\mathbf{x},t)& \\
	&\rho \mathbf{u}(\mathbf{x},t)=\sum_{i} \mathbf{c}_i f_i(\mathbf{x},t)&
	\end{eqnarray}
\end{itemize}


\section{Boundary conditions}
Boundary conditions are of great importance for obtaining any meaningful result. Most of the boundary conditions are given by macroscopic quantities like pressure, velocities, temperature and so on. Lattice Boltzmann Method is developed with concept of particle distribution functions in Mesoscopic scale. Therefore, this hydrodynamic quantities should be translated into conditions for particle distribution. This chapter will briefly introduce a few kind of boundary conditions that are widely used nowadays.

\subsection{Periodic Boundary Condition}
Periodic boundary condition could be the simplest boundary condition for LBM. In the periodic system, the fluid flow out through one face will reenter to the opposite face of the domain. So the edges of the simulation domain could be treated as they attach with the other side of itself. For the boundary nodes, their neighbouring nodes are located at the opposite side of the boundary. (Fig \ref{pbc})

\begin{figure}[!hpg]
\begin{center}
\input{pbc.pstex_t}
\end{center}
\caption{Periodic boundary condition}\label{pbc}
\end{figure}

Assume that the length of computational domain is $L$, periodic boundary condition could be describe in formula(take x direction as an example, the derivation for y direction is straightforward):
\begin{eqnarray}
&f_i(0,y,z,t+\Delta t)=f'_i(L,y,z,t)&\\
&f_i(L,y,z,t+\Delta t)=f'_i(0,y,z,t)&
\end{eqnarray}

\subsection{Bounceback Boundary condition}
Due to its simplicity, versatility and powerful capability to deal with extremely complex Boundaries.Bounceback condition plays a major role int he LBM simulation. This boundary condition is usually used at fluid-solid interfaces due to its corresponding to the no-slip condition. This boundary condition could be illustrate in Figure \ref{bounceback}. The densities moving toward to the solid are bounced back in to fluid domain with the incoming directions. In D2Q9 model, Bounceback condition can be described in term of equations as:

\begin{eqnarray}
f_2=f_4,&f_5=f_7,&f_6=f_8
\end{eqnarray}

The standard bounceback condition places the boundary on the lattice nodes, and it can keep the conservation of mass and momentum, but the accuracy is first order, where LBM is of second order. Inamuro investigated that the error produced by bounceback condition is sufficiently samll if the relaxation parameter $\tau$ is close enough to 2. (Inamuro et al. 1995)\cite{inamuro1995} The bounce back conditions can be used without any influence on the order of the LBM, if $\tau$ is chosen within a suitable range. Furthermore, the bounce bcak condition is the most dfficient one for arbitrary complex geometries (Breuer et al. 2000)\cite{breuer2000} A lot of researcher contributed for this ongoing discussion (Noble et al. 1995a,b; Ziegler, 1993; Skordos, 1993l Inamuro et al. 1995).A second order scheme is proposed by Ziegler (1993). Solid boundary is placed between two neighbouring lattice sites with same distance $\Delta x /2$. It is illustrate as Figure(\ref{bounceback}).   


\begin{figure}[!h]
\begin{center}
\input{bounceback.pstex_t}
\end{center}
\caption{Bounceback condition}\label{bounceback}
\end{figure}

\subsection{Flux Boundary}
Flux boundary conditions or Von Neumann Boundaries enable us constrain the flux at the boundaries. It was proposed by Zou and He (1997)\cite{zhouhe}. A velocity vector $\mathbf{u}_0=(u_0,v_0)^T$ is specified at the boundary, and the density and pressure will be computed.Due to the symmetry, we only derive the equations for north side boundary showed in Figure \ref{fluxbc}, the derivations for other cases are straightforward.\\


\begin{figure}[!h]
\begin{center}
\input{fluxboundary.pstex_t}
\end{center}
\caption{Flux condition}\label{fluxbc}
\end{figure}


After the streaming step, distribution function $f_0,f_1,f_2,f_3,f_5,f_6$ is known, and $f_4,f_7,f_8,\rho$ it to be determined. Suppose the velocity on the boundary is,

\begin{equation}
\mathbf{u}=\begin{pmatrix}
u_x \\ u_y
\end{pmatrix}
\label{boundaryv}
\end{equation}

The macroscopic density and velocity formula are:

\begin{eqnarray}
&\rho=\sum_{i} f_i& \\
&\mathbf{u} = \frac{1}{\rho}\sum_{i}f_i \mathbf{e}_i&
\label{macrorhov}
\end{eqnarray}

Formula \ref{macrorhov} and boundary velocity \ref{boundaryv} give three equations:

\begin{eqnarray}
&f_4+f_7+f_8=z\rho-f_0-f_1-f_2-f_3-f_5-f_6& \label{bc1}\\
&f_8-f_7=\rho u_x-(f_1-f_3-f_6+f_5) &  \label{bc2}\\
&f_4+f_7+f_8=\rho u_y+ (f_2+f_5+f_6) &\label{bc3}
\end{eqnarray}

We need another equation to solve all the four unknowns. It can be written by assuming that the bounceback condition holds in the direction normal to the boundary.

\begin{equation}
f_2-f^{eq}_2=f_4-f^{eq}_4
\label{zhouhebc}
\end{equation} 

Equation \ref{bc1} and \ref{bc3} lead to 
\begin{equation}
\rho=\frac{1}{1-u_y}[f_0+f_1+f_3+2(f_2+f_5+f_6)]
\label{bcdensity}
\end{equation}

We proceed by substituting equation \ref{bcdensity} into equation \ref{bc1},\ref{bc2} and \ref{bc3}, solve for $f_4,f_7,f_8$:

\begin{eqnarray}
&f_4=f_2+\frac{2}{3}\rho u_y& \\
&f_8=f_6-\frac{1}{2}(f_1-f_3)+\frac{1}{2}\rho u_x+\frac{1}{6}\rho u_y& \\
&f_7=f_5+\frac{1}{2}(f_1-f_3)-\frac{1}{2}\rho u_x+\frac{1}{6}\rho y_y&
\end{eqnarray} \\

Flux velocity boundary are widely used in the simulation with introduction of flux, it has special advantage in the application of multiphase flow simulation. It enable us to inject multiphase fluid with specified velocity simultaneity without any extra adjustment. Since the density of fluid will be calculate based on the number density of neighbouring nodes and boundary velocity. The phase status will be identified automatically. 

\subsection{Pressure/Density boundaries}
Pressure/Density boundaries or Dirichlet boundaries enable us to constrain the pressure/density on the boundaries during the simulation. Since there is an equation of state to relate pressure and density in LBM, constraining the pressure is equivalent to specify the density on the boundaries. Similar with the flux boundaries, after streaming, there are three unknown distributions plus one velocity, (We assume that the velocity tangent to the boundary $u_x$ is zero, and just compute the normal velocity $u_y$) the velocities will be computed based on specified density and distribution of neighbouring sites. The other unknowns will be obtained by solving the equation systems. The derivation of is very similar from above procedure, we just give the solution:

\begin{eqnarray}
&f_4=f_2-\frac{2}{3}\rho_0 u_y& \\
&f_7=f_5+\frac{1}{2}(f_1-f_3)-\frac{1}{6}\rho u_y& \\
&f_8=f_6-\frac{1}{2}(f_1-f_3)-\frac{1}{6}\rho u_y&
\end{eqnarray}


\subsection{Incorporating force term}\label{forcet}
The lattice Boltzmann equation given above is without the force term, however, in most applications, an external force such as gravity or interaction force between particles is applied to the fluid. Shan and Chen (1993)\cite{shanchen} proposed a form by changing the equilibrium velocity.Assume a foce $\mathbf{F}$ is applied. With Newton's law, write it as:

\begin{equation}
\mathbf{F}=m\mathbf{a}=m\frac{d\mathbf{u}}{dt}
\label{nl2}
\end{equation}

The mass is proportional to the density and relaxation time $\tau$ is the time of collision. We rewrite the equation \ref{nl2} to:

\begin{equation}
\Delta \mathbf{u}=\frac{\tau \mathbf{F}}{\rho}
\end{equation}

A new equilibrium velocity $\mathbf{u}^{eq}$ is used:

\begin{equation}
\mathbf{u}^{eq}=\mathbf{u}+\Delta \mathbf{u}=\mathbf{u}+\frac{\tau \mathbf{F}}{\rho}
\end{equation}

In the original paper of Shan and Chen, the macroscopic velocity was computed without any revision. In 1995, they presented a new paper pointing that, the macroscopic velocity need to be revised as the average of the velocity before the collision and the velocity after collision\cite{shanchen1995}:

\begin{equation}
\rho \mathbf{u}=\frac{1}{2}\rho (\mathbf{u}+\mathbf{u}')=\sum_{i} \mathbf{c}_if_i+\frac{\Delta t}{2}\mathbf{F}
\end{equation}


There are some other methods like adding force term in the evolution equations to match the corresponding macroscopic equations to the Navier-Stokes equation using the Chapman-Enskog expansion, or change the pressure term in the equilibrium distribution.A representation of the forcing term in shich discrete lattice effects are considered was proposed by Guo et al(2002)\cite{guo2002}.The scheme is based on D2Q9 lattice and the force is taken into account by adding an additional term to the LBE:

\begin{equation}
f_i(x+e_i\Delta t, t+\Delta t)-f_i(x,t)=-\frac{1}{\tau}[f_i(x,t)-f^{eq}_i(x,t)]+\Delta t F_i
\end{equation}

where the equilibrium equation $f^{eq}$ is defined as follow:

\begin{eqnarray}
&f^{eq}_i(\rho, u^*)=\omega_i \rho[1+\frac{e \cdot u^*}{c^2_s}+\frac{u^*u^*:(e_ie_i-c_s^2I)}{2c_s^4}]&\\
\text{with}&\rho u^*=\sum\limits_{i} e_if_i+mF\Delta t& \label{macroguo}
\end{eqnarray}

In equation (\ref{macroguo}) $m$ is a constant to be determined. The forcing term $F_i$ can be written in a power series in the particle velocity (Ladd and Verberg 2001)\cite{ladd2001}

\begin{equation}
F_i=\omega_i[A+\frac{B \cdot e_i}{c^2_s}+\frac{C:(e_ie_i-c^2_sI)}{2c^4_s}]
\end{equation}

Where $A,B$ and $C$ are functions of $F$ to be etermined by requiring that the moments of $F_i$ are consistent with the hydrodynamic equations. 
\begin{eqnarray}
&\sum\limits_{i} F_i=A& \\
&\sum\limits_{i} e_iF_i=B& \\
&\sum\limits_{i} e_ie_iF_i=c_s^2AI+\frac{1}{2}[C+C^T]& 
\end{eqnarray}

In order to study macroscopic behaviour, a multiscaling analysis using Chapman-Enskog expansion is used. Expanding $f_i{f+e_i\Delta t, t+\Delta t}$ about $x$ and $t$, and applying multiscaling expansions to the resulting continuous equations, compare the coefficient with Navier-Stokes equations,$A,B$ and $C$ can be determined and lead to the following force term:

\begin{equation}
F_i=(1-\frac{1}{2\tau} \omega_i [\frac{e_i-v}{c^2_s}+\frac{(e_i \cdot v)}{c^4_s}e_i]\cdot F)
\end{equation} 

Where $v$ is the fluid velocity defined by

\begin{equation}
\rho v=\sum\limits_{i} e_if_i+\frac{\Delta t}{2}F
\end{equation}

Kinetic viscosity is 

\begin{equation}
\nu=(\tau-\frac{1}{2})c^2_s \Delta t
\end{equation}

The macroscopic equations obtained by Chapman-Enskog expansion is given as
\begin{eqnarray}
&\frac{\partial \rho}{\partial t}+\nabla \cdot \rho v =0&\\
&\frac{\partial}{\partial t}(\rho v)+\nabla \cdot (\rho vv)=-\nabla p+ \nu \nabla \cdot [\rho (\nabla v+\nabla v^T)]+F&
\end{eqnarray}

The analysis above shows that, Guo force term can lead to exact Navier-Stokes equations.Guo and Zheng (2002)\cite{guo2002} analysis a few kinds of method in details and compare the simulation results.The force scheme proposed by Shan-Chen (1995)\cite{shanchen1995} is examined. If we neglect the term $O(\delta t^2)$ in momentum equation, the macroscopic equations corresponding to this method are

\begin{align}
\frac{\partial \rho}{\partial t}+\nabla \cdot \rho v =&0\\
\frac{\partial}{\partial t}(\rho v)+\nabla \cdot (\rho vv)=&-\nabla p+ \nu \nabla \cdot [\rho (\nabla v+\nabla v^T)]+F\\
&\nabla \cdot[\rho \nu^2(\nabla a+\nabla a^T)]
\end{align}

Since the the term $\nabla \cdot[\rho \nu^2(\nabla a+\nabla a^T)]$ is not negligible, this method is mainly used for flows exposed to a constant body force.

Guo and Zheng (2002)\cite{guo2002} analysis a few kinds of method in details and compare the simulation results.


\section{Multiphase flow}
LBM has been proved that it is a promising technique in simulating multiphase problem. That is the true strength of LBM. Especially for the problem with complex phase interactions with complex geometries. And it is capable for various phenomenons like evaporation, condensation, cavitation, capillarity that are difficulty for traditional CFD technique.Since the phase separation due to the long-range interaction between the moldcules of the fluid, therefore, additional term should be added to the Naavier-Stokes equations for traditional CFD technique.While the interface should also the treated separately. The traditional CFD technique using VOF, level set to capture interface evolution, but they are only capable to describe interface movement in big scale for few cases, and can not capture the interface with small size. Due to the nature of LBM, interactions between different phases or component could be described with ease. This feature enable us to find the key to preserve the essential physic of multiphase flow.\\


The first LBM multiphase flow model is proposed by Gunstensen et al (1991)\cite{color1991} with the name color function model. Two components represent two types of fluid with their own distribution functions, and follow their own LBM equation. They are named as red particles and blue particles. Collision step include self-interactions and cross-interactions with other types of particles. Colour function gradient was introduced to calculate the surface tension between different phases. To segregate the phases, mixing near the interface should be minimised, a procedure called recoloring is proposed for the minimising process. This procedure is a very time costing step, and this model also have some numerical stability problem for high density ratio and large surface tension value.\\


The free energy model was developed by Swift et al (1996) \cite{freeenergy}. This model includes the thermodynamic equilibrium of phases, and a term describing the surface tension is added to the equilibrium distribution function. This character make free energy model easier to specify the surface tension value than other multiphase models.However, lack of Galilean invariance is a principle shortcoming of free energy model.\\

\subsection{Shan-Chen model}
We will work with Shan and Chen model (1993), which is a pseudo potential model developed by Shan and Chen. The principal characteristic of this model is an interaction force between particles is introduced to have a consistent treatment of equation of state of non-ideal gas. The interactions are involved with nearest neighbours, so the model is computationally efficient.For D2Q9 model,an attractive force between particles is described as:

\begin{equation}
F(\mathbf{x},t)=-G\Psi (\mathbf{x},t)\sum_{i=1}^{8}w_i \Psi (\mathbf{x}+\mathbf{e}_i\Delta t,t)\mathbf{e}_i
\label{Shanchenforce}
\end{equation}

$G$ is a parameter determine the interaction strength between neighbouring particles. It also determine if the interaction is attraction or compulsion. In multi-component flow, the same component particles will attract each other while particles from different component will repulse each other. The component separation can be obtained. And $\Psi$ is the interaction potential with the form:

\begin{equation}
\Psi (\rho)=\Psi_0exp(\frac{-\rho_0}{\rho})
\label{potentialfunction}
\end{equation}

Where $\Psi_0$ and $\rho_0$ are constant parameters. In Shan-Chen model, interaction forces between molecules are
accounted, and the interaction forces change the original ideal gas
equation of state into a non-ideal equation of state. With a
non-ideal equation of state, phase separation will take place when
the temprerature falls below the critical temperature of the fluid.
For multi component fulid simulation, intre-species interations are
also considered. \\

Expanding equation \ref{Shanchenforce} at $x$, the interaction force can be written as:
\begin{equation}
F=-c^2_sG\Psi \nabla \Psi -\frac{1}{2}\Delta x^2c^2_sG\Psi \nabla \nabla^2 \Psi +o(\Delta x^3)
\label{expansionforce}
\end{equation}

Chapman-Enskog expansion shows the relation between lattice Boltzmann equation and Navier Stokes equation, and the force term can be combined with the pressure gradient term. To include the interaction force term, the pressure should be modified  and results in a non-ideal EOS as:

\begin{equation}
p(\rho)=\rho c^2_s +\frac{1}{2}Gc^2_s \Psi ^2(\rho)
\label{sukopeos}
\end{equation}


This interaction force term is applied using method mentioned in Chapter 3. Application of this force term leads to a non-ideal gas equation of state:

\begin{equation}
P=\rho RT+\frac{GRT}{2}[\Psi(\rho)]^2
\end{equation}

Where $RT=\frac{1}{3}$ is used in EOS.\\


The second term is related with interface stress corresponding to
the gradient free energy. It control the geometry of interface as
well as the surface tension. It is dependent on the mesh resolution
due to the existing of $\Delta x^2$.

In thermodynamic theory, the equilibrium pressures and densities
with particular temperature could be obtained by using Maxwell
construction (Rowlinson and Widom 1982). The Maxwell Construction
can be stated as
\begin{equation}
\int _{V_{m,l}}^{V_{m,v}} PdV_m = P{V_{m,v}-V_{m,l}}
\end{equation}

\begin{figure}[!h]
\begin{center}
\includegraphics[width=3in]{maxwellconstruction.eps}
\end{center}
\caption{The molar volumes of the liquid and vapor and vapor pressure can be determined by Maxwell Construction}\label{MxC}
\end{figure}



The area under the curve must equal the area of the rectangle
defined by vapor and molar volumes and the vapor pressure.Fig (\ref{MxC})\\

The Applying the mechanical equilibrium condition on the stress
tensor across the interface determine the interface structure.

\begin{equation}
\sigma= \int_{-\infty} ^{\infty} (P_{zz}-P_{xx})dz
\end{equation}

Assumption of interface lying in the x-y plane is made for above
equation.And $P_{zz}$ and $P_{xx}$ are the stresses normal and
parallel to the interface respectively. The value of surface tension
in the interaction potential  model have been reported by
Shan(2008)\cite{shan2008}:

\begin{equation}
\sigma=-\frac{eGc_s^4 e_4 \Delta x^2}{2} \int_{-\infty} ^{\infty}
(\frac{d\Psi}{dz})^2 dz
\end{equation}

From above analysis, we can find that, the surface tension is
related with the chosen of potential function which is equivalence
to the EOS, and as well as the lattice geometry. In order to make
the surface tension independent of  of EOS

By introducing a new term to instead of the second term in equation
(\ref{expansionforce}) It's possible to make the surface tension independent
from the selection of EOS. Change the weight of the term $\psi
\nabla \nabla ^2 \psi$ with a coefficient $\kappa$

\begin{equation}
F^s = \kappa \psi \nabla \nabla ^2 \psi
\end{equation}

The discretization form of above function could be obtained by using
the pair-wise interaction between lattice site:

\begin{equation}
F^s (x) = \kappa \psi(x) \sum _{y} \lambda_{xy} \psi(y) |x-y| /
\Delta x^3
\label{addShanChenterm}
\end{equation}

A larger set of lattice sites in addition to the nearest neighbors
are need to obtain higher order difference operator for the
potential function $\psi$. The coefficient $\lambda_{xy}$ can be
determined by Taylor expansion of function
(\ref{addShanChenterm}).As a result, the total interaction force can
be written as:

\begin{equation}
F=F^0 + F^s \approx -\nabla (\frac{1}{2}c_s^2 G \Psi ^2) + \kappa
\Psi \nabla \nabla ^2 \Psi
\end{equation}
Shan-Chen pseudo potential model could be improved also for fluid-surface interaction. In this case, the force should be described is interaction force between the fluid particles and solid particles. Neighbouring node potential function will not be computed, instead, a parameter $G_{ads}$ is introduced to give the interaction strength between the fluid and solid, all the neighbouring nodes that are solid need to be account for.

\begin{equation} 
F_{solid}(\mathbf{x},t)=-G_{ads} \Psi (\mathbf{x},t)\sum_{i} w_i s(\mathbf{x}+\mathbf{e}_i\Delta t,t)\mathbf{e}_i
\end{equation}

Where function $s$ indicate if the node is a solid node: it takes value 1 if the node $(\mathbf{x}+\mathbf{e}_i\Delta t)$ is solid, and 0 if the node is fluid. The applying of the force is similar as fluid interaction force term. The long-range interaction does not affect the shear viscosity (Shan and Doolean 1995)\cite{shanchen1995}


\subsection{Spurious current and higher order isotropy multiphase LBE model}
The spurious velocity is the unphysical circulating velocity near
the interface. The presence will affect the accuracy of the
simulation and also cause numerical stability problems. This
phenomenon is a common problem for many other multiphase simulations
like VOF , color gradient LBM and free energy LBM. The magnitude is
found related with the width of the interface and the surface
tension. Increased surface tension will lead to big spurious
velocity while wider interface gives smaller spurious velocity. The
existence of spurious limited ability of the pseudo potential model
to simulate high density ratios or low viscosities. Shan (2008)
reported that the spurious is due to the lack of sufficient isotropy
in the calculation of the gradient term for force. Finite difference gradient of force with higher degree of isotropy
is also introduced.\\

Long-range intermolecular interaction can be introduced by defining
a interaction potential between neighboring site:

\begin{equation}
V(x,x')=G(x,x')\psi (x) \psi (x')
\end{equation}

Where $\psi=\psi(\rho)$ is the "effective mass" that defines the
details of the interparticle interaction. The force can be obtained
by

\begin{equation}
F=-G(|e_a|)\psi(x) \sum\limits_a \psi (x+e_a)e_a
\label{forcesc}
\end{equation}

The momentum can be showed conserved due to the summation over all neighbours interacting with $x$. The right hand side of equation (\ref{forcesc}) is the finite difference representation of $-\psi \nabla \psi$. The interaction strength can be determined by more neighbouring sites with different distance. The calculation of gradient of potential function could achieve sufficient degree of isotropy and accuracy. Analysis of finite difference scheme could help us determine the value of coefficients $G(|e_a|)$ (Shan 2006)\cite{shan2006}. Expand $\psi(x+e)$ as Taylor series on $x$

\begin{equation}
\psi(x+e)=\sum\limits_{n=0}^{\infty} \frac{1}{n!}[\nabla ^{(n)}\psi (x)]\cdot \underbrace{ee\ldots e}_{n}
\end{equation}

The approximation of the gradient of potential function is obtained:
\begin{align}
\sum\limits_{a=1}^{d} w(|e_a|^2)\psi (x+e_a) e_a=&(\nabla \psi)\cdot E^{(2)}+\frac{1}{3!}(\nabla^{(3)} \psi) \cdot E^{(4)} \\
&\frac{1}{5!}(\nabla ^{(5)} \psi) \cdot E^{(6)}
\end{align}

where the weight $E^{(n)}$ is defined as
\begin{equation}
E^{(n)}_{i_1i_2\cdots i_n}=\sum\limits_{a} w(|e_a|^2)(e_a)_{i_1} \cdots (e_a)_{i_n}
\end{equation}

By choosing appropriate value of weight$w(|e_a|)$, unit $E^{(2)}$ and isotropic $E^{(n)}$ can be obtained. The solution for $w(|e_a|)$ which yield isotropic $E^{(n)}$ is given by Shan (2006)\cite{shan2006}


\begin{table}[htbp]
\begin{center}
\begin{tabular}{lcccccccc}
\toprule

& Tensor & w(1) & w(2) & w(3) & w(4) & w(5) & w(6) & w(8) \\
\midrule


&$E^{(4)}$ & 1/3 & 1/12 & & & & & \\
2D &$E^{(6)}$ & 4/15 & 1/10 & &1/120 & & & \\
&$E^{(8)}$ & 4/21 & 4/45 & & 1/60 &2/315 & & 1/5040 \\
\\
&$E^{(4)}$ & 1/6 & 1/12 & & & & & \\
3D &$E^{(6)}$ & 2/15 & 1/15 &1/60 &1/120 & & & \\
&$E^{(8)}$ & 4/45 & 1/21 &  2/105 &5/504 &1/315 & 1/630 & 1/5040 \\
\bottomrule

\end{tabular}
\end{center}
\caption{Weights that yield unit $E^{(2)}$ and isotropic $^{(n)}$}
\end{table}


\subsection{Discussion}
Due to the its clear physical background of modeling fluid by considering the collision, propagation, and interactions between fluid particles, LBM is a promising method in simulation of multiphase flow, multicomponent flow with various thermodynamical phenomenons which is difficult for traditional CFD technique like capillarity, condensation, evaporation, surface tension and cavitation. The description of particle dynamics enable LBM to deal with extremely complex boundaries like porous media. On the other hand, LBM is also an efficient, it avoid solving the pressure Poisson equation which is the most time-consuming step in traditional CFD. Instead the equation of state is used in LBM, which is better to describe the physics of multiphase flow. And as only the local information is used in the computing, LBM is intrinsically suitable for parallel computation. \\


In spite of huge number of reasons of choosing this promising methods, there are still a lot of restrictions for its application. As LBM is developed from continues Boltzmann equations which is used for compressible fluid, the Navier-Stokes equations can be recovered when pressure variation is small and low Reynolds number. In addition, numerical instability problem will emerge when the viscosity is small. Although the LBM has an extraordinary strength for multiphase and multicomponent flow simulation, the density ratio between different phases is constrained below around 50. High density ratio will cause numerical stability problem. Improvement can be obtained by using appropriate equation of state for each component and modified interaction force term between different component. The numerical stability and limit of viscosity can be improved by using Multi-Relaxation time LBM. Higher order of isotropy of gradient operator can reduce the spurious current as well as improve the numerical stability.

\subsection{Kinetic Model for Multiphase flow}
The continuous Boltzmann equation has mainly been used to solve
supersonic flows (K,Xu 1994)\cite{xu1994} This is due to the extreme complexity
of the collision kernel when dealing with dense, interacting
particles and the high requirement for computer resources. The
original Shan-Chen potential model only include attraction force
between particles, however, when the gas is compressed, the
repulsive force will dominate the gas model. Enskog theory take into
account the effects of molecular volume on molecular transport
properties.However, the long-range molecular interaction among
molecules is ceglected in Enskog's theory, we use mean-field theory
to describe the long-range molecular interactions. (He, Shan and
Doolen 1998) proposed a kinetic model that combines Enskog's theory
for dense fluids and mean-field theory for long-range molecular
interaction.

Continuum Boltzmann equation with BGK collision model is used to
built the kinetic model:

\begin{equation}
\frac{\partial f}{\partial t}+ \xi \cdot \nabla f+ F \cdot
\nabla_{\xi}f = -\frac{f-f^{eq}}{\lambda}
\end{equation}

Where $f(x,\xi,t)$ is the distribution function in the phase space
$(x, \xi)$, $\xi$ is the microscopic velocity and $F$ is the an
external body force. $f^{eq}$ is the Maxwell-Boltzmann distribution
function.

The macroscopic quantities are written as:

\begin{eqnarray}
&\rho=\int f d \xi& \\
&\rho \mathbf{u}= \int \xi fd\xi&\\
&\frac{D}{2} \rho RT= \int (\xi- \mathbf{u})^2 f d\xi&
\end{eqnarray}

The term $\nabla _{\xi} f$ cannot be computed directly since the
distribution function is dependent on the microscopic velocity which
is unknown. An approximation of using equilibrium distribution
function is made to compute this term:

\begin{equation}
\nabla _{\xi}f \approx \nabla_{\xi} f^{eq} = -\frac{\xi -
\mathbf{u}}{RT}f_{eq}
\end{equation}

Consequently, simplified Boltzmann equation is obtained as:

\begin{equation}
\frac{\partial f}{\partial t}+ \xi \cdot \nabla f=
-\frac{f-f^{eq}}{\lambda}+\frac{F \cdot (\xi -u)}{RT}f^{eq}
\label{equationenskog}
\end{equation}

There are two important factors for nonideal gas: the intermolecular
attraction and the exclusion volume of molecules. In He Shan and
Doolen model, the intermolecular attraction is treated using
mean-field approximation. Averaged force potential is used to
describe the intermolecular attraction force:

\begin{equation}
V(r_1)=\int _{r_{12}> \sigma} u_{attr}(r_{12})\rho(r_2)d r_2
\end{equation}

Where $r_{12}= |r_1-r_2|$ and $u_{attr}(r_{12})$ is the attractive
component of the intermolecular pairwise potential. $\sigma$ is the
diameter of the molecules. Expanding $\rho(r_2)$ about $r_1$ and
with the assumption that $\nabla \rho$ is small, approximation of
$V$ is achieved with first two leading term:

\begin{equation}
V=-2a\rho -\kappa \nabla ^2 \rho
\end{equation}

Where $a$ and $\kappa $ is constant given by:

\begin{eqnarray}
&a=-\frac{1}{2}\int _{r>\sigma} u_{attr} (r) d r&\\
&\kappa =\frac{1}{6} \int_{r>\sigma} r^2 u_{attr} (r) d r&
\end{eqnarray}

The effect of the exclusion volume of the molecules coule be
approximated by(S.Chapman and T.G.Cowling, the mathematical theory
of non-uniform gases 1970):

\begin{equation}
-f^{eq}bp\chi (\xi-u) \cdot \nabla ln(\rho^2 \chi)
\label{exclusionv}
\end{equation}

where $b=2\pi \sigma^3/3m$ $m$ is the mass of a single particle.
$\chi$ is the increase in collision probability due to the increase
in fluid density:

\begin{equation}
\chi = 1+\frac{5}{8}b\rho + 0.2869(b\rho)^2 +0.1103 (b\rho)^3+
\ldots
\end{equation}

Combine the equation (\ref{exclusionv}) with the intermolecular
attraction term, a equivalent force term is obtained:

\begin{equation}
F=-\nabla V-b\rho RT\chi \nabla ln(\rho^2 \chi)
\end{equation}

Chapman-Enskog analysis yield the corresponding macroscopic momentum
equation as:

\begin{equation}
\frac{\partial (\rho u)}{\partial t}+ \nabla \cdot (\rho uu)=-
\nabla \cdot P' +\nabla \cdot \tau + \rho a
\end{equation}

Where $\tau$ is the viscous stress $\tau=\mu_E S$ viscosity $\mu_E$
is defined as:

\begin{equation}
\mu_E = \rho c_s^2 (\frac{\tau}{\chi}-\frac{1}{2})\delta_t
\end{equation}

$P'$ is the thermodynamic pressure tensor with the definition as:
\begin{equation}
P'_{\alpha \beta} = p\delta_{\alpha \beta}+\kappa \frac{\partial
\rho}{\partial x_{\alpha}}\frac{\partial \rho}{\partial x_{\beta}}
\end{equation}
\begin{equation}
p=\rho RT(1+b\rho \chi)-a \rho^2 -\kappa \rho \nabla^2 \rho
-\frac{\kappa}{2} |\nabla \rho|^2
\end{equation}

The term $\nabla \cdot P'$ can be written as:
\begin{equation}
-\nabla \cdot P'=-\nabla p_0 +\kappa \rho \nabla \nabla^2 \rho
\end{equation}

The first term $\nabla p_0$ on the right hand side is the gradient
of equation of state:

\begin{equation}
p_0=\rho RT(1+b\rho \chi)-a \rho^2
\end{equation}

The second term $\kappa \rho \nabla \nabla^2 \rho$ represent the
surface tension. Different equation of state can be employed by
modifying the term value $p_0$.The force term can be written as:

\begin{equation}
F_t=\rho a_t = -\nabla \psi +F_s+\rho a \label{forcetermenskog}
\end{equation}

$\psi$ is a function related with equation of state:
\begin{equation}
\psi=p_0-p_{id}=p_0-\rho RT
\end{equation}


To solve equation \ref{equationenskog} numerically, we first
discrete the equation in time:
\begin{equation}
f(x+\xi \delta t,\xi, t+\delta t)-f(x,\xi,t)=-\int _t^{t+\delta t}
\Lambda (f-f^{eq}) dt+\int _t^{t+\delta t} \frac{F \cdot
(\xi-u)}{RT}f^{eq} dt
\end{equation}

The first integration is assumed to be constant over one time step,
and a trapezoidal rule is used to compute the integrand of the
second term

\begin{equation}
f(x+\xi \delta t,\xi, t+\delta t)-f(x,\xi,t)=-\Lambda (f-f^{eq})
\Delta t+ \frac{1}{2}(\frac{F \cdot (\xi -u)}{RT}f^{eq}|_{t}+\frac{F
\cdot (\xi -u)}{RT}f^{eq}|_{t+\Delta t})
\end{equation}




The LBE in the general form is:

\begin{equation}
f_i(x_i+c_i\Delta t, t+\Delta t)-f_i(x_i,t)= -\sum_j
\Lambda_{ij}(f_j-f_j^{eq})|_{(x,t)}+\frac{1}{2}[S_i
(x,t)+S_i(x+c_i\Delta t, t+\Delta t)]\label{MRTenskog}
\end{equation}

where $c_i$ is the discrete velocity, $S_i$ represent the force term
in Eq. (\ref{forcetermenskog}) The equation (\ref{MRTenskog}) has a
implicit term which usually requires inverse calculation. To
maintain the explicit scheme, we need eliminate this implicity, an
new variable is introduced to do so

\begin{equation}
\bar{f}_i =f_i-\frac{1}{2}S_i \Delta t \label{newtranformation}
\end{equation}

With the transformation in equation (\ref{newtranformation}),
equation (\ref{MRTenskog}) can be converted into an explicit scheme
as:

\begin{equation}
\bar{f}_i(x+c_i\Delta t, t+\Delta t)-\bar{f}_i(x,t)= -\sum_j
\Lambda_{ij} (\bar{f}_j-f^{eq}_j)|_{(x,t)}+\sum_j
(I_{ij}-\frac{1}{2}\Lambda_{ij})S_j|_{(x,t)}\Delta t
\end{equation}

Where $I$ is an identity matrix. The kinematic viscosity $\nu$ is
defined as
\begin{equation}
\nu=(\frac{1}{s_8}-\frac{1}{2})c_s ^2 \Delta t
\end{equation}

And the bulk viscosity is related with $s_2$:
\begin{equation}
\zeta =(\frac{1}{s_2}-\frac{1}{2})c_s^2 \Delta t
\end{equation}
 The macroscopic quantities can be calculated from the summations:
\begin{eqnarray}
&\rho=\sum\limits_{\alpha} f_{\alpha}=\sum\limits_{\alpha} \bar{f}_{\alpha}&\\
&\rho u=\sum\limits_{\alpha} f_{\alpha}e_{\alpha}=\sum\limits_{\alpha}
\bar{f}_{\alpha}e_{\alpha}+\frac{1}{2}F\delta t&
\end{eqnarray}




\section{Multi-Relaxation time LBM}
LBGK use single relaxation time and a Taylor series expansion of the
Maxwell-Boltzmann equilibrium distribution function for the
equilibria. There are some limitations for viscosity, low viscosity
will cause numerical instability. There is no conclusion for the
reason of this numerical instability, some have attributed them to
the nonexistence of an H theorem. (B.Boghosianl 2001)\cite{boghosianl2001}.
There exist two school of thought on how to overcome numerical
instabilities(McCracken 2005)\cite{mccracken2005} The first suggests using nonpolynomial
equilibrium distribution functions, however, doing this will
significantly increase the computing requirements and will not
always lead to correct Naiver Stokes equations. The others propose
the use of multiple relaxation times in order to reduce the
numerical instabilities. \\

MRT allows independent adjustment of bulk and shear viscosities which enable Multi-Relaxation time LBM significantly improve the numerical stability for low viscosity fluid. For Shan-Chen LBGK model, the numerical stability problem partly from the spurious velocity across the interface caused by the decrease of viscosity. The collision term is relaxed by single parameter $\tau$ in LBGK model, while, it could be instead by a matrix $\Lambda$:

\begin{equation}
f_i(x+c_i \Delta t)-f_i(x,t)=-\Lambda_{ij}[f_j-f_j^{eq}], \quad i=1,2,\ldots,b
\end{equation}

Matrix $\Lambda$ is a full matrix with constant. LBGK model could be obtained by specifying $\Lambda$ as a diagonal matrix with identical value:

\begin{equation}
\Lambda_{ij}=\frac{1}{\tau}\delta_{ij}
\end{equation}

The macroscopic quantities are calculated in the same way as LBGK model.Instead of treating distribution function, MRT employs several moment corresponding to macroscopic quantities and their flux. These quantities can relax with different time scales. A matrix $M$ transform the distribution functions $f_i$ from distribution space into moment space. 

\begin{equation}
m=M\cdot f, \quad f=M^{-1}\cdot m
\end{equation}

Collision is carried out in moment space,multiply the transformation matrix $M$, the left and right hand side of equation \ref{collisionmrt} can be transformed into the moment space:

\begin{equation}
f_i'(x,t)=f_i(x,t)-\Lambda_{ij}[f_j-f_j^{eq}]
\label{collisionmrt}
\end{equation}
\begin{equation}
\mathbf{m}'=\mathbf{m}-\mathbf{S}[\mathbf{m}-\mathbf{m}^{eq}]
\end{equation}

Where $\mathbf{m}^{eq}=\mathbf{M}\mathbf{f}^{eq}$ is the equilibrium equation in moment space. $S=M \Lambda M^{-1}=diag(s_1,s_2,\ldots,s_b)$. The corresponding relax time for moment $m_i$ is $s^{-1}_i$. After collision step, the moment $m'$ is transformed into distribution function space by multiplying $M^{-1}$ for streaming step which will be carried out in the same way as LBGK model.


For D2Q9 model, the transform matrix $M$ is given as \cite{lallemand2000}:
\begin{equation}
M=\begin{pmatrix}
1  &1  &1  &1  &1  &1  &1  &1  &1  \\
-4 &-1 &-1 &-1 &-1 &2  &2  &2  &2\\
4  &-2 &-2 &-2 &-2 &1  &1  &1  &1\\
0  &1  &0  &-1 &0  &1  &-1 &-1 &1\\
0  &-2 &0  &2  &0  &1  &-1 &-1 &1 \\
0  &0  &1  &0  &-1 &1  &1  &-1 &-1\\
0  &0  &-2 &0  &2  &1  &1  &-1 &-1\\
0  &1  &-1 &1  &-1 &0  &0  &0  &0\\
0  &0  &0  &0  &0  &1  &-1 &1 &-1
\end{pmatrix}
\end{equation}

and the corresponding moments are:
\begin{equation}
m=(\rho,e,e^2,j_x,q_x,j_y,q_y,p_{xx},p_{xy})^{T}
\end{equation}

Where $e$ is the energy, $j_x,j_y$ are the momentum in $x$ and $y$ direction, $q_x,q_y$ are energy flux, $p_{xx},p_{xy}$ is the stress tensor.\\


The model is athermal, therefore only density and momentum is conserved. Collision do not change the conserved quantities. Energy is not conserved in athermal model. Since collision has no effect on conserved quantities, the relaxation time will be 0 for these quantities, and an assuming that the nonconserved moments relax linearly towards to their equilibrium that are functions of conserved quantities. The relaxation parameters and equilibrium functions of moments are:



\begin{equation}
S=(0,s_e,s_{e^2},0,s_q,0,s_q,s_{\nu},s_{\nu})
\end{equation}
\begin{equation}
m^{eq}=\rho(1,-2+3u^2,\alpha+\beta u^2,u_x,-u_x,u_y,-u_y,u^2_x-u^2_y,u_xu_y)^{T}
\end{equation}

Where $\alpha$ and $\beta$ are parameters, and kinematic viscosity and volume viscosity are:
\begin{eqnarray}
&\nu=c^2_s(\frac{1}{s_{\nu}}-\frac{1}{2})\Delta_t& \\
&\zeta=c^2_s(\frac{1}{s_{e}}-\frac{1}{2})\Delta_t&
\end{eqnarray}





\section{Numerial simulation result}
\subsection{Phase separation}
With Shan-Chen pseudo potential model, we can simulate phase separation and its dynamics.A $200 \times 200$ domain with average density of $200 mu lu^{-1}$ is computed. The initial density is specified with random variations to prevent a metastable situation.Parameter $\Phi _0=4$ and $\rho_0=200$ and $G=-120$ are used for the calculation of interaction force term.The surface area volume is minimized due to the consequence of free energy minimization. In the liquid-vapor system, condensation and evaporation phenomenons can also discovered during the simulation. This can be due to the relatively high vapor density.

\begin{figure}[h]
\begin{center}
\includegraphics[width=5in]{phaseseperation.eps}
\end{center}
\caption{Liquid-vapor phase separation dynamics at time 0, 100, 200,400,2000,6000 $ts$}
\end{figure}

\subsection{Maxwell Construction}
Flat interfaces are exceptionally important because the vapor pressure above interface at equilibrium is the saturation vapor pressure. The pressure across the flat interface should be zero according to the Laplace law. Density of liquid and vapor yield from the simulation can be compared with Maxwell Construction result. A 50*200 domain was set to compute the yield density from Shan-Chen model. EOS in equation (\ref{sukopeos}) with $\Psi_0=4, \quad \rho_0=200$ and $G=-120$ is computed. Shan-Chen multiphase model is used to separate the phase, two force incorporating scheme was used, one is the shifting velocity scheme proposed by Shan and Chen (1995) \cite{shanchen1995}, it has been widely used in LBM multiphase simulation, another force term is presented by Guo (2002)\cite{guo2002} with consideration of discrete lattice effects. Simulation results are compared with numerical solution from Maxwell Construction. 

\begin{table}[!htbp]
\begin{center}
\begin{tabular}{lccc}
\toprule

& $\rho_{liquid}$ & $\rho_{vapor}$ & $|u_{max}|$  \\
\midrule
Shan-Chen & 524.4 & 85.75 & 1.6291E-4 \\
Guo force term & 514.7 &79.69 & 1.1204E-4 \\
Analytical solution &514.64 & 79.705 & \\
\bottomrule

\end{tabular}
\end{center}
\caption{Liquid and vapor density yield from Shan-Chen model with two different force terms and the numerical solution from Maxwell Constructions. All values are in lattice unit}
\end{table}

The results show that liquid and vapor density yield from Shan-Chen model with original force term has a significent difference from analytical solution, while Guo force term can lead to very accurate result. The way of incorporating force term also have big effects on spurious velocity $u_max$. Guo force term also gives smaller spurious velocity across the interface. The analysis in Chapter \ref{forcet} shows Shan-Chen original force term is a first order scheme, which can only recover mass equation, and has a error term $\nabla \cdot[\rho \nu^2(\nabla a+\nabla a^T)]$ in momentum equation. Guo's force term can recover both mass and momentum equation. It may be the reason of leading to more accurate liquid and vapor density compared with Shan-Chen force term.
  

\subsection{Surface tension}
By simulating a series of drops of different size and measuring the inside and outside pressure, the surface tension can be estimated with Laplace law. The slop of 1/radius vs. pressure difference will be the surface tension. The following numerical examples using improved Shan-Chen model (\ref{expansionforce1a}) which enable the surface tension independently from the EOS. 

\begin{equation}
F=-c^2_sG\Psi \nabla \Psi -\frac{1}{2}\Delta x^2c^2_sG\Psi \nabla \nabla^2 \Psi +o(\Delta x^3)
\label{expansionforce1a}
\end{equation}

A non-ideal EOS(\ref{sukopeos2}) with the potential function defined in Equation(\ref{potentialfunction}) is used. 
\begin{equation}
p(\rho)=\rho c^2_s +\frac{1}{2}Gc^2_s \Psi ^2(\rho)
\label{sukopeos2}
\end{equation}

Pressure field of different size drops and different $\kappa$ are compute.  

\begin{figure}[!htbp]
\begin{center}
\includegraphics[width=3in]{surface_tension.eps}
\end{center}
\caption{Surface tension estimation with different $\kappa$}
\label{sur1}
\end{figure}


\begin{figure}[!htbp]
\begin{center}
\includegraphics[width=3in]{surface_tension2.eps}
\end{center}
\caption{Surface tension $\sigma$ change along with parameter $\kappa$}
\label{sur2}
\end{figure}

The Fig (\ref{sur1}) shows the Laplace law is satisfied with $\kappa=12,20,30$. And different $\kappa$ can result different surface tension $\sigma$ which is calculated with Laplace law. Surface tension resulted from small surface tension can fit the Laplace law better, while result from big surface tension produce bigger error. Fig (\ref{sur2}) shows the linear relation between paramter $\kappa$ and surface tension $\sigma$. This numerical example demonstrate the validation of the new Shan-Chen model with improvement which enable us to control the surface tension apart from equation of state.





\subsection{Spurious current}
The relation between spurious current and degree of isotropy is examed by simulating a single stationary drop in vapor with different order of isotropic gradient operator. Numerical example from Shan (2006)\cite{shan2006} showed that higher order of isotropic gradient operator can significantly reduce the spurious velocity near the interface. Since Shan-Chen model with original force term will give relatively big yield density compared with analytical solution, numerical simulation is made to study if the use of higher order isotropic gradient operator can reduce this unexpected density difference.\\

\begin{figure}[!htbp]
\begin{center}
\includegraphics[width=3in]{spv2.eps}
\end{center}
\caption{Reduction of spurious velocity with different order isotropy gradient tensor. Higher order isotropy can significantly reduce the magnitude of spurious velocity}
\label{sur2}
\end{figure}


The effect of using higher order isotropy gradient tensor for the equilibrium density is also studied. The result below shows the equilibrium density of flat interface test.




\begin{table}[!htbp]
\begin{center}
\begin{tabular}{ccc}
\toprule

Order of isotropy & $\rho_{liquid}$ & $\rho_{vapor}$   \\
\midrule
4 & 524.44  & 85.734 \\
6 & 524.892 & 86.0416 \\
8 & 525.735& 86.6422 \\
Analytical solution & 514.64 & 79.705  \\
\bottomrule

\end{tabular}
\end{center}
\caption{Equilibrium density with different degree of isotropy.Inclusion of high order isotropy force term can not improve the resulting equilibrium density.}
\end{table}




The improved Shan-Chen model, or the Enskog dens gas model need to compute the term $\nabla \psi$ separately. Since the surface tension term $\rho \nabla \nabla^2 \rho$ is normally small compare with pressure term, the discrete implementation of this term will affect the numerical stability of the simulation. Several different order of accuracy finite difference of the term  $\nabla \psi$ are computed for the simulation of a single stationary drop in vapor, the spurious velocities near interface are compared.

\begin{figure}[!htbp]
\begin{center}
\includegraphics[width=3in]{spv1.eps}
\end{center}
\caption{Reduction of spurious velocity with different order accuracy of approximation of $\nabla \psi$. The magnitude of spurious velocity is decreasing significantly due to introduction of the higher order accuracy term}
\label{sur2}
\end{figure}


The equilibrium density yield from model using different approximation scheme is studied.Second, forth, and sixth order accuracy of term $\nabla \psi$ is used in the simulation.


\begin{table}[!htbp]
\begin{center}
\begin{tabular}{ccc}
\toprule

Order of accuracy & $\rho_{liquid}$ & $\rho_{vapor}$   \\
\midrule
2  & 524.888 & 86.0172 \\
4  &523.538 & 85.1135 \\
6 & 523.282 & 84.9454 \\
Analytical solution & 514.64 & 79.705  \\
\bottomrule

\end{tabular}
\end{center}
\caption{Equilibrium density with different accuracy order scheme of $\nabla \psi$. No significant improvement is obtained in this test. The effect of order of accuracy for equilibrium density is limited.}
\end{table}

\section{Literature review}
Shan and Doolen(1995) \cite{shanchen1995} revise the force term of original Shan-Chen model. Macroscopic velocity will be obtained by average the velocity before the collision and after. Detailed analysis for Shan-Chen model was carried out, mutual diffusivity can be obtained as a function of concentrations of the components. It does not depend on the fluid velocity and is fully Galilean invariant. However, this model suffered from the spurious velocity at the interface. \\

Several efforts were made to improve the stability of LBM for multiphase simulation. Inamuro et al(2004)\cite{inamuro2004} simulate two-phase flows with large density difference by solving the pressure Poisson equation with free energy based LBM. Density ratio of 1000 was obtained. Solving pressure Poisson equation is very time-consuming and lose the thermodynamic meaning of the problem. Yuan and Schaefer (2006)\cite{yuan2006} compare several types of equation of state and investigate the effects of interaction potentials on the density ratio range in single component simulation. By applying appropriate interaction potential which related with equation of state, density ratio of 1000 can be obtained. Their study showed that Peng-Robinson(P-R) and Carnahan-Starling (C-S) EOS can reach density ratio over 1000. The improvement of numerical stability comes from the reduce of sound speed. However, their study is restricted to single component model.\\


A forcing term in which discrete lattice effects are considered is proposed by Guo and Zheng(2002)\cite{guo2002}. Chapman-Enskog analysis shows that exact Navier-Stokes equations can be obtained.Several other existing force treatments are also studied and compared with this force scheme.It is found that none of these method match the Navier-Stokes equations in the general case.


Multiple-relaxation-time (MRT) lattice Boltzmann was developed by d'Humieres (1992) \cite{dhumieres1992}. The main idea of MRT is use several relaxation time parameters for different moments of macroscopic quantities. And the collision step is carried out in moment space, while streaming step is still done in velocity space as LBGK. MRT model can overcome several defects of the LBGK model such as fixed Prandtl number (Pr=1 for the BGK model) and fixed ratio between the kinematic and bulk viscosities. And the stability improvement by using the MRT scheme would reduce the computational effort by at least one order of magnitude while maintaining the accuracy of the simulations.(d'Humieres et al 2001)\cite{dhumieres2001}. The numerical stability of MRT is studied by (Lallemand and Luo 2000)\cite{lallemand2000} in detailed.Their analysis showed the MRT obtained a big improvement of numerical stability compared with LBGK model. MRT model for incompressible flow is developed by Du et al (2006) \cite{du2006}. Rectangular grid based multiple-relaxation-time lattice Boltzmann equation was developed by Bouzidi et al(2001)\cite{bouzidi2001}\\

Morton number $Mo=\frac{g\rho^3 \nu^4}{\sigma^3}$ is introduced to show the ratio of viscosity and surface tension. \\


Due to the interaction are limited with attraction force only, when
the gas is compressed to near its 'hard sphere' volume, the
repulsive force will dominate the van der Waals gas model.
Thermodynamic theorem is satisfied only when the potential function
is in the form of exponential function:

\begin{equation}
\Psi (\rho)=\psi e^{-\rho_0/ \rho}
\end{equation}

Therefore the original Shan-Chen model has very limited ability to
simulate any desired EOS. Incorporating with Enskog dens gas theorem
can resolve EOS limitation. Several LBM incorporate with Enskog dens
gas theorem has been developed. (He and Doolen 2002, Martys 2001,
Luo 2000, Luo 1998)\cite{he2002} \cite{martys2001} \cite{luo2000} \cite{luo1998}. \\

A discrete model based on the Boltzmann equation with a body force was developed by He,Shan and Doolen(1998)\cite{he1998}. The Boltzmann equation is discretized in a way that preserves the derivation of the hydrodynamic equation from the Boltzmann equation. A mean-field approximation is used in their model for the interparticle interaction.They reported that, in this type of model, anisotropy might be inevitable if a nonlocal interaction is not included in momentum space. \\

LBM multiphase models based on Enskog theorem have been extended with MRT algorithm in 2D and 3D(McCracken,2005, Premnath, 2008)\cite{mccracken2005} \cite{premnath2008}. Improvement in stability and capability for simulating low viscosity fluid was reported to be achieved.\\

Comparison of different multiphase LBM model in thermodynamic aspect has been given by He and Doolen (2002)\cite{he2002}. Detailed analysis of thermodynamic foundations of kinetic theory was given. Demonstration shows that kinetic model is consistent with thermodynamic theory and can lead to correct not only mass and momentum equations but also energy transfer equation.\\

The spurious problem was first studied by Shan (2006) \cite{shan2006}. Analysis show that the spurious current near interface is due to the insufficient isotropy of the discrete gradient operator. Finite difference gradient operators with higher order of isotropy was proposed and the spurious current was found decreasing significantly.  Sbragaglia (2007)\cite{sbragaglia2007} analysis the flows near interface in details and extended pseudopotential method is developed with using multirange pseudopotential. The extended pseudopotential model also permits to tune the equation of state and surface tension independently of each other.\\









\bibliographystyle{ieeetr}
\bibliography{ref}
\end{document} 

